function y = tridiag(c,b,a,f)
% function z = tridiag(c,b,a,f)
%   this function is a simple triadiagonal solver
%   for Ay = f, where A is a tridiagonal matrix
%   consisting of upper, lower, and main diagonals
%   of vectors c, b, and a, respectively.
N = length(a);
w = zeros(N,1);
v = w; z = w; y = w;

% set w1 = ad1
w(1) = a(1);
v(1) = c(1)/w(1);
z(1) = f(1)/w(1);

c = [c' 0]';  b = [0 b'];

for i = 2:N
    w(i) = a(i)-b(i)*v(i-1);
    v(i) = c(i)/w(i);
    z(i) = (f(i)-b(i)*z(i-1))/w(i);
end

y(end) = z(end);

for i = (N-1):-1:1
    y(i) = z(i)-v(i)*y(i+1);
end


